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^ I We present a straightforward integration method to compute the abundance and tem- 

00 ' perature evolution in explosive scenarios. In this approach the thermal equation is im- 

^ ■ plicitely coupled with chemical equations in order to avoid instabilities and ensure a gentle 

Q . transition from the normal combustion regime to the quasi (QSE) and complete nuclear 
statistical equilibrium (NSE). Two nuclear networks, with 14 nuclei (a-network) and 86 

^ I nuclei (including protons and neutrons) respectively, have been considered. The scheme 

^ I is suitable to cope with a variety of explosive burning regimes. 

O ' 1. Introduction 

^ I Violent phenomenae in astrophysics are often related with the sudden release of nuclear 
binding energy. In some cases, when material is very degenerate, such release is practically 
isochoric and we refer it as explosive ignition. In other situations there is some expansion 
and the evolution takes place at nearly constant pressure and high temperature. In both 
5^ I combustion regimes, however, the nuclear system evolves from the initial composition to 
a final state characterized by the coexistence of hundreds or even thousands elements in 
chemical equilibrium. From the point of view of stellar modelization, different burning 
regimes translates into a variety of subroutines in order to properly handle all combustion 
stages. It is frequent to use a second order implicit scheme to integrate the chemical 
equations and assume NSE wherever density and temperature are high enough. While 
this approach has proven useful in many situations, it would be worthwile an integration 
algorithm suitable to handle rapid combustion, including the limiting cases of QSE and 
NSE, fast and stable enough to keep large time-steps. A possible solution is to couple 
molar fractions with temperature through the energy or entropy equation as proposed 
by E.Miiller [Q. Taking a small network of 13 nuclei, from a to ^^Ni, he showed how 
the inclusion of temperature effects was decisive to avoid the numerical oscillations which 
otherwise show up during and after the relaxation to the NSE regime. In this paper we 
extend the previous work of Miiller considering a larger nuclear network and providing 
some more test cases. In addition we also propose a slightly different integration method 
which is very robust and straightforward to implement. Nowadays, the effort to write 
the proposed numerical scheme for large networks (even 'kilo-nuclide' networks) can take 
great advantage of libraries of theoretical reaction rates recently published @ . 
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2. Description of the method 

First of all the chemical equations are written in their implicit form and linearized 
taking +^ = + AF^ as described in [|]. In addition we also take T"+i = + AT 
and make a first order Taylor expansion of nuclear rates with respect the temperature. For 
instance, for photodesintegration rates we take A"^"^ = Aj(T")+(?Ai/c?T|j'nAT. Neglecting 
second order terms like AFj AF^ or AF^AT the discretized system of equations reads: 



^ = Y. ruYr^Yk - ( E n.YJ" I AF, + ^ A^ AF„ - A, AF, + ^ r^F," AF^ 

k,l \ j I m k,l 

- Y r^.Y-AY, + [Y ^k^kY^ - E + E 4^;: - A:F," ] AT 

+ E rkm'' - E r^fY-yr + E >^mYZ - KYl^ (1) 

k,l j rn 

where = pNa{a, v)ij stands for particle reactions, Aj is for photodesintegrations, an 
upper comma means to take the temperature derivative and the remaining symbols have 
the usual meaning. In addition we need an extra equation to close the system. We took 
the energy equation, which for an adiabatic process can be written as follows: 



E5^.AF,-|(^ 



Eventual diffusion terms should be added to the independent terms to the right of 
Equation 2. This method can be applied to arbitrarly large networks. The computational 
effort to integrate the system with and without thermal coupling is practically the same. 

2.1. Tests realized 

The nuclear evolution of dense plasmas under extreme thermal conditions involve a large 
quantity of nuclei. However, in current multidimensional simulations of supernovae and 
others explosive events the size of the network is strongly constrained by computational 
efficiency. Small networks of around ten nuclei are used in this case, which approximately 
give the correct nuclear energy rate. Thus, we have checked the method with two networks 
of different size: an a-chain of 14 nuclei from a to ^^Zn (including ^^C, and ^^O binary 
reactions) and a network of moderate size with 86 nuclei from protons and neutrons to 
^^Zn. Nuclear rates were taken from and the sparse matrix solver was that of 0]. 

The following cases of astrophysical interest were successfully calculated: a) isochoric 
combustion of an initial fuel composed by O at p = lOl^g.cm'^ followed by an 

adiabatic expansion in the hydrodynamical time (Figure 1, upper right), b) photodesin- 
tegration of a system composed of ^^Fe at constant density of lQ'^g.cm~^ (bottom left) 
and, c) isobaric propagation of a nuclear flame in a medium composed by ^^C +^^0 and 
p = 10^ g.cm~^ (bottom right). As a comparation, the isochoric burning case was also 
calculated without thermal coupling. In this case the oscillations of both, abundances and 
temperature, are clearly visible when NSE is approached (upper left, using 14 nuclei). 
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Figure 1. Nuclear combustion for the four cases discussed in the text. 



Fortunately, once the thermal coupling was included the QSE and the NSE regimes 
were smoothly achieved with no fluctuation in all cases, allowing the scheme to take 
longer time-steps in these phases. More details about the performance of the method will 
be presented elsewhere 0. 
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